load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

;#HERE YOU WILL PUT THE PATH TO control (a2h2) and lowice (a2g4) ensemble members

filsw1= systemfunc ("ls "+"/Users/mypath/a2g4/ua/ua_allmon_levs_djf*.nc")
filsc1= systemfunc ("ls "+"/Users/mypath/a2h2/ua/ua_allmon_levs_djf*.nc")
filsw2= systemfunc ("ls "+"/Users/mypath/a2g4/va/va_allmon_levs_djf*.nc")
filsc2= systemfunc ("ls "+"/Users/mypath/a2h2/va/va_allmon_levs_djf*.nc")

fw1    = addfiles (filsw1, "r")
fc1    = addfiles (filsc1, "r")
fw2    = addfiles (filsw2, "r")
fc2    = addfiles (filsc2, "r")

ListSetType (fw1, "join")
ListSetType (fc1, "join")
ListSetType (fw2, "join")
ListSetType (fc2, "join")

print(fw1)
print(fc1)
print(fw2)
print(fc2)

;*******************************************
;Read in variables
;*******************************************
lat2d=fc1[0]->latitude
lon2d=fc1[0]->longitude

lev=2 
ens =10
steps=60
uempcss =  new((/ens*steps,256,512/),float)
vempcss =  new((/ens*steps,256,512/),float)
uempc =  new((/ens,steps,256,512/),float)
vempc =  new((/ens,steps,256,512/),float)


do i=0,ens-1

time=fc1[i] ->time
ntime= dimsizes(time)
print("ntime   =  "+ntime)
uempsc =  new((/ntime,3,256,512/),float)
uempsc = fc1[i]->ua
vempsc =  new((/ntime,3,256,512/),float)
vempsc = fc2[i]->va

last=ntime-1
first=ntime-60
uempc(i,:,:,:)=uempsc(first:last,lev,:,:)
vempc(i,:,:,:)=vempsc(first:last,lev,:,:)

do m=0,steps-1
im=(i*steps)+m
uempcss(im,:,:)=uempc(i,m,:,:)
vempcss(im,:,:)=vempc(i,m,:,:)
end do

delete(time)
delete(ntime)
delete(last)
delete(first)
delete(uempsc)
delete(vempsc)
end do

uempwss =  new((/ens*steps,256,512/),float)
vempwss =  new((/ens*steps,256,512/),float)
uempw =  new((/ens,steps,256,512/),float)
vempw =  new((/ens,steps,256,512/),float)

do i=0,ens-1

time=fw1[i] ->time
ntime= dimsizes(time)

uempsc =  new((/ntime,3,256,512/),float)
uempsc = fw1[i]->ua
vempsc =  new((/ntime,3,256,512/),float)
vempsc = fw2[i]->va

last=ntime-1
first=ntime-60
uempw(i,:,:,:)=uempsc(first:last,lev,:,:)
vempw(i,:,:,:)=vempsc(first:last,lev,:,:)

do m=0,steps-1
im=(i*steps)+m
uempwss(im,:,:)=uempw(i,m,:,:)
vempwss(im,:,:)=vempw(i,m,:,:)
end do

delete(time)
delete(ntime)
delete(last)
delete(first)
delete(uempsc)
delete(vempsc)
end do

uempcss!0 = "alltimesteps"
uempwss!0 = "alltimesteps"

vempcss!0 = "alltimesteps"
vempwss!0 = "alltimesteps"

printVarSummary(uempcss)
printVarSummary(uempwss)

;************************************************
; streamfunction derivation
;************************************************

 sfunc_c=uv2sfvpF(uempcss,vempcss)
 sfunc_w=uv2sfvpF(uempwss,vempwss)

 delete(uempcss)
 delete(vempcss)
 delete(uempwss)
 delete(vempwss)

 sfc=sfunc_c(0,:,:,:)
 sfw=sfunc_w(0,:,:,:)

 sfc!0 = "alltimesteps"
 sfw!0 = "alltimesteps"
 sfc!2 = "x"
 sfw!2 = "x"
 sfc!1 = "y"
 sfw!1 = "y"

 printVarSummary(sfc)
 delete(sfunc_c)
 delete(sfunc_w)
;************************************************
; average over the simulations/years
;************************************************

  finfo_tcs = sfc(0,:,:)
  finfo_tcs = dim_avg(sfc(y|:,x|:,alltimesteps|:))     ; reorder but do it only once [temporary]
  finfo_tws = sfw(0,:,:)
  finfo_tws = dim_avg(sfw(y|:,x|:,alltimesteps|:))

;************************
; calculate the anomaly:
;************************

  finfo_ano=finfo_tws
  finfo_ano=(finfo_tws -finfo_tcs)/100000.

  finfo_ano!0 = "lat2d"
  finfo_ano@lat2d = fc1[0]->latitude

  finfo_ano!1 = "lon2d"
  finfo_ano@lon2d =fc1[0]->longitude

;________________
;write data

name1 ="SF250hPa_control_djf_ensmean"
outf1 = addfile(name1+".nc","c")
outf1->time = time
outf1->lat  = lat2d 
outf1->SF250   = finfo_tcs

name2 ="SF250hPa_lowice_djf_ensmean"
outf2 = addfile(name2+".nc","c")
outf2->time = time
outf2->lat  = lat2d
outf2->SF250  = finfo_tws

;-------------------------
; making the plots
;-------------------------
wks = gsn_open_wks("pdf","SFano_250hPa_djf_ecearth3p")

gsn_define_colormap(wks,"BlueDarkOrange18")

plot = new(1,graphic)                          ; create a plot array

  res          = True
  res@gsnDraw  = False                          ; don't draw
  res@gsnFrame = False                          ; don't advance frame, these 2 are going to be set with panel plot - gsn_panel!
  res@cnInfoLabelOn = False                     ; turn off cn info label
  res@cnLinesOn           = True 
  res@cnFillOn            = True          ; turn on color
  res@lbLabelBarOn        = True           ; turn off individual cb's
 
  res@cnLevelSelectionMode =  "ExplicitLevels"

 res@cnLineLabelsOn=False
 res@mpCenterLonF     = 290. 

res@gsnPolar   = "NH"
res@mpMinLatF            = 0.

res@cnMonoLineColor = True 
res@cnMonoLineThickness= True
res@cnLineThicknessF = 2.
res@cnLineColorF= "black"

res@gsnContourNegLineDashPattern=1
res@gsnContourZeroLineThicknessF=3

res@mpGeophysicalLineThicknessF = 2.

res@gsnLeftString = "DJF"
res@gsnLeftString = ""
res@gsnCenterString = ""

 res@cnFillColors = (/2,3,4,5,8,10,0,0,13,14,15,16,17,18/)
 
 res@lbAutoManage = False
 res@lbLabelFontHeightF  = 0.017
 res@pmLabelBarHeightF  = 0.08

res@cnLevels =(/-20.,-15.,-10.,-7.5,-5.,-2.5,0,2.5,5.,7.5,10.,15.,20./)

res@gsnRightString = ""
plot(0) = gsn_csm_contour_map(wks,finfo_ano(:,:),res)

;________________________
; creating the panel:
;________________________
 resP            = True

resP@gsnPanelLabelBar    = False;True                ; add common colorbar
resP@lbLabelFontHeightF  = 0.007               ; make labels smaller

resP@gsnPanelBottom   = 0.05                   ; add space at bottom
;resP@gsnPanelFigureStrings= (/"a)","b)","c)"/) ; add strings to panel
res@txFontHeightF     = .24


resP@gsnMaximize    = True                ; maximize plot
gsn_panel(wks,plot,(/1,1/),resP)


